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1  Introduction 

The  problem  addressed  in  this  project  is  that  of  tracking  a  single  satellite  in  the  neighborhood  of  several 
other  closely  spaced,  similar  satellites  where  there  is  measurement  mixing  between  satellites  observed  at  a 
ground  based  sensor.  This  situation  has  been  referred  to  as  the  “furball  problem.” 

Tracking  a  single  object  or  a  group  of  objects  in  formation  are  both  well  studied  and  understood  processes. 
What  makes  the  “furball  problem”  unique  is  that  the  problem  of  tracking  a  single  object  in  a  formation  is  a 
generally  undeveloped  capability,  especially  with  the  type  of  sensor  considered  for  this  particular  case.  The 
ground  based  sensor  considered  takes  a  single  measurement  per  sampling  period.  The  single  measurement 
can  potentially  be  from  any  of  the  satellites  in  the  sensor’s  measurement  region.  This  implies  that  the  sensor 
returns  a  single  stream  of  measurements  that  are  “mixed”  in  the  sense  that  the  satellite  from  which  the 
measurement  originated  is  not  known  to  the  sensor. 

The  purpose  of  this  report  is  to  formally  address  the  ‘  furball  problem”  and  present  several  effective  data 
association  methods  for  finding  a  viable  solution.  As  will  be  shown,  these  data  association  techniques  can 
be  used  to  effectively  sort  out  measurements  and  provide  reliable  tracking. 

The  report  is  structured  in  the  following  manner.  The  problem  is  formally  defined  and  the  various  data 
association  techniques  used  to  attack  this  problem  arc  briefly  presented  in  Sections  2  and  3.  Section  4 
provides  a  comparison  of  the  computational  complexity  between  existing  methods  of  tracking  through  batch 
processes  and  iterative  methods  described  in  Section  3.  In  Section  5  we  provide  an  in  depth  description 
of  the  types  of  orbits  being  considered  along  with  a  description  of  the  sensor  based  frames  from  which  the 
measurements  are  taken.  Simulation  results  are  presented  in  Section  6,  comparing  the  performances  of  the 
algorithms  discussed  in  Section  3.  Finally,  some  conclusions  and  possible  areas  of  future  work  are  discussed 
in  Section  6. 


2  Problem  Definition 

The  exact  problem  considered  in  this  report  is  a  system  consisting  of  j  satellites  in  either  the  same  orbit  or 
in  orbits  very  close  to  each  other.  For  our  purposes,  two  orbits  being  “close”  can  be  defined  as  keeping  all 
the  orbital  parameters  the  same  for  each  orbit  while  perturbing  the  eccentricity  of  one  by  a  small  amount. 

A  ground  based  sensor  is  being  used  in  an  attempt  to  track  a  single  one  of  the  j  closely  spaced  satellites. 

The  model  used  to  describe  the  dynamics  of  the  satellites  written  down  in  an  earth  centered  inertial  (ECI) 
coordinate  frame  is 

£'i( k)  =  f[k,Xi(k  -  1)]  +Vi(k),i  =  1,2 (1) 

where  /  contains  the  dynamics  of  the  system,  and  the  additive  process  noise  Vi(k)  is  zero-mean  and  white 
with  covariance  Q(k).  Since  maneuvering  is  not  considered  in  this  initial  work,  the  dynamics  contain  no 
control  vector.  The  measurement  due  to  each  satellite  is 

Zi(k)  —  h[k,Xi(k  -  1)]  +Wi(k),i  ==  1,2,.  ..J  (2) 

where  the  additive  noise  Wi(k)  is  also  assumed  to  be  zero-mean  and  white  with  covariance  R(k ),  and  the  func¬ 
tion  h  contains  the  coordinate  transformation  between  ECI,  South/East/Z  (SEZ),  and  Range/ Azimuth/ Elevation 


(R/Az/El)  coordinates.  The  ECI  coordinates  are  Cartesian  coordinates  sueh  that  the  x-axis  is  aligned  with 
the  vernal  equinox,  the  2-axis  is  straight  up  through  the  north  pole,  and  the  y- axis  is  90  degrees  away  from 
the  x-axis  in  the  equatorial  plane  sueh  that  the  coordinates  are  right-handed.  The  SEZ  coordinates  are 
Cartesian  coordinates  based  at  the  latitude  and  longitude  position  of  whatever  sensor  is  being  used  to  take 
measurements,  the  x-axis  aligned  with  due  South,  the  y-axis  aligned  with  the  East,  and  the  2-axis  pointing 
direetly  out  tangent  to  the  surfaee  of  the  earth.  The  range,  azimuth,  and  elevation  are  measured  from  the 
SEZ  frame.  The  range  is  a  radial  position  measurement  from  the  sensor,  the  azimuth  is  the  angular  measure 
from  true  north  (— x  in  SEZ),  and  the  elevation  is  the  angular  measure  sueh  that  the  tangent  plane  at  the 
sensor  has  measure  zero  and  direetly  above  (2  in  SEZ)  has  measure  90  degrees. 

The  measurements  for  this  system  are  taken  in  R/Az/El  coordinates  and  a  single  measurement  is  taken 
at  eaeh  time  step  where  the  origin  of  that  measurement  is  not  known.  For  example,  consider  two  satellites  in 
elosely  spaced  orbits,  where  we  are  trying  to  traek  satellite  1.  Five  measurements  are  taken,  one  of  whieh  is 
randomly  from  satellite  2.  One  possibility  for  how  the  list  of  measurements  recorded  by  the  sensor  appears 
is 

(z(/c),  z(k  +  1),  z(k  +  2),  z(k  +  3),  z(k  +  4)}  =  { z\{k ),  z\{k  +  l),zi(fc  +  2),  z2(fc  +  3),  zi(k  +  4)}. 

Only  the  set  of  measurements  (with  unknown  origin)  is  provided.  The  aetual  identity  of  individual  measure¬ 
ments  on  the  right  hand  side  is  not  known.  The  data  association  algorithm  must  sort  out  whieh  measurement 
in  the  set  originated  from  the  satellite  of  interest. 

The  above  example  ean  be  used  to  introduce  uniformity  index  (UI),  a  parameter  in  this  study.  The 
uniformity  index  is  the  percentage  of  the  total  measurements  that  are  assumed  to  be  from  the  eorreet 
satellite  of  interest.  In  the  example,  one  of  the  five  measurements  is  always  assumed  to  be  wrong,  though 
its  position  remains  unknown.  Thus,  the  UI  for  this  ease  is  80%. 

Due  to  the  faet  that  the  sensor  is  considered  to  be  on  the  surface  of  the  earth,  in  reality  there  will  elearly 
be  times  when  the  two  satellites  moving  in  elose  proximity  will  be  “out  of  view.”  This  problem  will  be 
touehed  upon  in  the  future  work  seetion,  but  for  most  of  this  report  will  be  ignored.  Thus  for  this  report,  it 
is  assumed  that  there  is  only  a  single  sensor  being  used  to  traek  satellites  through  an  entire  orbit. 


3  Algorithms 

3.1  The  Batch  Filter 

The  Bateh  Filter  is  a  widely  used  filtering  technique  that  is  often  applied  to  systems  were  there  is  no  mixing 
of  measurements  and  where  estimation  does  not  need  to  be  done  in  an  online  manner.  The  Batch  Filter  is 
used  in  these  instances  beeause  it  will  in  general  converge  to  zero  RMS  error.  Seetion  A  on  computational 
complexity  will  illustrate  a  reason  the  Bateh  Filter  is  not  always  preferable. 

The  implementation  of  the  Bateh  Filter  is  as  follows.  Measurements  from  a  single  satellite  are  taken. 
The  continuous  time  dynamics  for  the  system  are  assumed  to  be  known  perfectly  and  ean  be  written 

x(t)  =  f{t,x(t)).  (3) 

Then  for  k  measurements  taken  at  discrete  time  steps  the  continuous  dynamics  are  used  to  solve  baekwards 
in  time  k  time  steps  to  obtain  k  estimates  of  the  state  at  time  zero,  i.e.,  (xi(0),  x2(0),  ...,Xfc(0))  where <Xi(0) 
is  obtained  by  integrating  baekwards  from  the  first  measurement  by  dt,  x2(0)  by  integrating  backwards  from 
the  seeond  measurement  by  2  •  d£,  and  so  on  up  to  k.  These  k  estimates  of  the  initial  state  are  then  averaged 
to  obtain  a  “best  estimate,”  x(0),  of  the  state  at  t  =  0.  The  dynamics  (3)  ean  again  be  used,  this  time  to 
integrate  forward  in  time,  starting  from  the  x(0)  and  integrating  forward  in  time  to  time  k.  This  value  is 
then  the  current  estimate  of  state.  This  process  is  repeated  for  time  k  +  1  up  to  the  final  measurement.  In 
principle,  if  the  noise  associated  with  the  measurements  is  white  with  zero  mean,  the  bateh  estimates  should 
eventually  converge  to  zero  RMS  error  as  k  goes  to  infinity. 

t 

3.2  The  Extended  Kalman  Filter 

The  EKF  is  a  well  known  filtering  method  for  nonlinear  systems  [  ].  The  EKF  is  an  estimation  algorithm 
that  maps  the  current  estimate  of  the  state,  x(k  —  1| k  —  1),  forward  one  time  step.  The  eurrent  estimate  is 
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conditioned  on  the  current  state  given  the  current  measurement  [  ] 

x(k  -  1 1 A;  -  1)  4  E[x(k  -  1  )\Zk~l] 

where 

zk~' ±{z(j),j  =  1, . ,k-  1} 

is  the  cumulative  set  of  measurements  up  to  time  k  —  1.  The  associated  state  error  covariance  matrix  is 

P(k  -  1| k  -  1)  =  E{[x(k  -  1)  -  x(k  -  1| k  -  l)][x(k  -  1)  -  x(k  -  1|A:  -  l)]'|Zfc_1}.  (4) 

We  proceed  with  finding  x(k\k  —  1),  the  prediction,  by  expanding  the  nonlinear  function  (1).  The 
expansion  is  accomplished  by  evaluating  the  Taylor  series  of  (1)  about  the  current  estimate:  x(k  —  1 1 A;  —  1 ) 


Tv 

x(k)  =  /[&, x(k  —  1| k  —  1)]  +  fx(k  ~  1)[^(A:  —  1)  —  x (k  —  l\k  —  1)]  +  1/2  ei[x(k  —  1)  —  x(k  —  1| k  —  1)]' 

i=  1 

•  fxx(k  —  l)[x(k  —  1)  —  x(k  —  1 1 A;  —  1)]  H-  higher-order  terms  +  v(k  —  1)  (5) 

where  e*  is  the  ith  Cartesian  basis  vector  and  fx(k  —  1)  is  the  Jacobian  of  the  vector  / 

fx(k  —  1)  =  [ SJxf  {k  ~  1)  X)]x=x(fc-l|/c-l)* 

Similarly,  the  Hessian  of  the  ith  component  of  /  is 

flxx(k  -  1)  =  [yx  y*  f(k  -  1,  *^)]x=x(/c— ij/c— i.)  • 

The  prediction  of  the  state  at  time  k  is  then  obtained  by  taking  the  expectation  of  the  Taylor  series  expansion 
(5)  conditioned  on  Zk~l  and  neglecting  higher-order  terms: 


72^ 

x{k\k  -  1)  =  f[k,  x(k  -  l\k  -  1)]  +  1/2  £  e^Tr  [fxx(k  ~  1  )P(k  -  1|  k  -  1)]. 

2  =  1 

The  state  prediction  error  is 


TLX 

x(k\k  —  1)  =  fx(k  —  1  )x(k  —  1|  k  —  1)  +  l/2^^ei[x'(k  —  1 1  A:  —  \)flxx(k  —  1  )x(k  —  1|  k  —  1) 

2=1 

-  WL(k  -  1  )P(k  -  1| k  -  1)]]  +  v(k  -  1) 

and  its  covariance  is 


nx  nx 

P(k\k  -  1)  4  E[x(k\k  ~  l)x'(k\k  -  l)\Zk~ J]  =  fx(k  -  1  )P(k  -  1| k  -  1  )f'x  +  1/2  ££  e^TrifUk  -  1) 

2  =  1  J  =  1 

*  P(k  -  II k  -  1  )fxx(k  -  1  )P(k  -  1| k  -  1)]  +  Q(k  -  1). 

The  measurement  prediction  is 

nz 

z(k\k  —  1)  =  h[k,  x(k\k  —  1)]  +  1/2  ejTi[hxx(k)P(k\k  —  1)]. 

2=1 

The  innovation  is 


v(k)  =  z{k)  —  z{k\k  —  1) 

=  h[k,  x(k\k  —  1)]  T  w(k), 


3 


and  the  innovation’s  covariance  is 


S(k)  4  hx(k)P(k\k  -  1  )h'x(k)  +  lfrjrjreie'jlhUQPiklk  -  1)  •  hix(k)P(k\k  -  1)]  +  R(k) 

2=1  j— 1 

where 

hx(k)  4 

and 

Kx{k)  =  [Vx  Vx  hi(k,x)}'x=mk-iy 

The  EKF  gain  is 

W{k)  4  p(k\k  -  1  )tix{k)S~l{k). 

The  filtered  estimate  is  then 

x(k\k)  =  x(k\k  -  1)  +  W{k)v{k)  (6) 

and  the  filtered  estimate  error  covariance  is 

P(k\k)  =  P(k\k  -  1)  -  P(k\k  -  1  )tix(k)S~l(k)hx(k)P(k\k  -  1) 

=  [/-  W(k)hx(k)]P(k\k-l). 


3.3  Gating 

Gating  is  a  method  of  determining  whether  or  not  it  is  likely  that  a  particular  measurement  came  from  the 
object  being  tracked.  Since  we  are  assuming  that  the  measurement  at  each  time  step  is  normally  distributed 
around  the  truth  measurement,  it  is  possible  to  define  a  region  in  the  measurement  space  where  there  is  a 
high  probability  of  finding  the  measurement 

U( 7)  =  {z  ■■  u'(k)S~l(k)v{k)  <  7}.  (7) 


The  probability  of  gating  the  truth  measurement  is 

f7  yi./2-l  exp|_7/2} 


-f 


2nzl2T(nz/2) 


CZ7, 7  >  0 


where  T(-)  is  the  gamma  function  and  nz  is  the  dimension  of  the  measurement  space.  Because  the  mea¬ 
surement  error  is  assumed  Gaussian,  this  norm  is  chi-square  distributed.  If  the  measurement  error  is  not 
Gaussian,  such  as  in  the  case  of  accounting  for  the  underlying  geometry  of  the  system  in  the  measurement, 
gating  can  still  be  carried  out  with  some  modification  [  ]. 

As  defined  in  Sections  1  and  2,  at  each  time  step  there  is  exactly  one  measurement  of  unknown  origin  that 
is  received.  If  the  measurement  falls  within  the  gate,  it  is  passed  on  for  use  in  the  EKF.  If  the  measurement 
falls  outside  the  gate,  no  measurement  is  passed  to  the  EKF,  and  the  EKF  simply  predicts  forward  one  time 
step. 


3.4  Probabilistic  Data  Association 

Probabilistic  Data  Association  has  been  used  when  multiple  measurements  are  received  at  each  time  step 
in  the  tracking  of  an  object  in  a  region  with  uniformly  distributed  environmental  clutter  [  ].  Given  this 
assumption  (uniformly  distributed  clutter),  the  probability  of  each  measurement  Zi  being  the  true  mea¬ 
surement  from  the  object  of  interest  is  computed.  The  PDA  is  used  in  conjunction  with  a  filtering  algorithm, 
such  as  the  EKF,  where  updated  estimates  Xi(k\k)  due  to  each  measurement  Zi  are  computed  and  the  overall 
updated  state  estimate  is 


mk  , 

x(k\k)  =  ]xj(k\k)Pi(k ),  (8) 

2  =  0 

where  each  i*  is  of  the  form  (G),  rrik  is  the  number  of  measurements  at  time  k ,  and 
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^A(fc)  =  1, 

i=0 

i.e.,  the  probabilities  are  mutually  exclusive  as  well  as  exhaustive.  When  tracking  in  cluttered  environments 
where  multiple  measurements  are  received  at  each  time  step,  gating  is  often  used  to  limit  the  computational 
complexity.  When  gating  is  used,  rrik  in  (8)  is  the  number  of  gated  measurements.  Thus  for  the  system 
considered  in  this  report,  rrik  —  1  always. 

When  tracking  in  cluttered  environments  with  some  types  of  sensors  (e.g.,  radar  or  sonar),  the  probability 
Pd  of  detecting  the  object  of  interest  is  often  less  than  1.  That  is,  at  each  time  step,  the  probability  Pd 
that  one  of  the  multiple  measurements  received  is  due  to  the  object  of  interest  is  less  than  1.  For  this  report, 
Pd  =  UI. 

The  i  =  0  case  in  (8)  represents  the  hypothesis  that  no  true  measurement  from  the  object  of  interest  is 
received  at  time  k.  This  implies  that  for  the  type  of  sensor  described  in  Section  2,  if  the  one  measurement 
per  time  step  is  determined  to  not  be  from  the  object  of  interest,  the  update  (8)  will  be  entirely  based  on 
the  i  =  0  case. 

The  error  covariance  associated  with  (8)  is 

P(k\k)  =  0o(k)P(k\k  -  1)  +  [1  -  P0{k)}Pc{k\k)  +  P(k) 


where  P(k\k  —  1)  is  the  same  as  is  in  the  EKF, 


where 


P(k)  =  W(k) 


mk 

PJ0i(k)i'i(k)i''i(k)  -  u(k)u'{k) 

_l—  1 


W'(k) 


Ui{k)  =  Zi(k)  -  z(k\k  -  1), 

rnk 

=  'P&i(k)vi(k), 

2—1 


and 

Pc(k\k)  =  [I-  W(k)H(k)\P(k\k  -  1). 


Clutter  measurements  are  assumed  to  be  uniformly  distributed  throughout  the  tracking  volume  with 
density 

\  _  mk 
Vk' 

where  Vk  is  the  volume  of  the  validation  region 


=  CnMk)  |1/2  =  Cnzln‘,2\S(k)\l/2 


and  nz  is  the  dimension  of  the  measurement  and  cUz  is  the  volume  of  the  n^-dimensional  unit  hypersphere. 
The  f3i(k )  probabilities  are  computed  as  [ 


where 


and 


0i(k) 


b  +  ET= iV 


Po  (k) 


b  +  ET^ej 


i  -  1 


=  exp 


i 'i{k)S  l{k)vi{k) 

2 


b  ^  A|27r5(fc)|1/2(1  -  PdPg)/Pd 
=  (27r/7r/2Al4/Cn,(l  -  PdPg)/Pd. 
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3.5  Kolmogorov-Smirnov  Tests 

Kolmogorov- Smirnov  (KS)  tests  provide  another  approach  for  determining  the  probability  of  whether  a  mea¬ 
surement  in  a  set  originated  from  the  object  of  interest  [  ].  Consider  a?i,  £2, xn,  independent  observations 
of  a  random  variable  with  unknown  cumulative  distribution  function  (CDF)  F(x).  If  the  hypothesis  is 

H  :  F(x)  =  F0(x), 

then  any  test  of  this  hypothesis  is  a  goodness-of-fit  test  [  ].  This  hypothesis  states  that  the  theoretical  CDF 
of  the  random  variable  is  equal  to  the  CDF  from  which  actual  measurements  are  drawn.  The  KS  tests  and 
many  simple  variants  are  goodness  of  fit  tests.  Empirical  CDFs  are  formed  for  a  window  of  n  measurements 
and  some  metric  is  then  applied  to  measure  the  distance  between  theoretical  and  empirical  CDFs. 

For  n  ordered  samples 

£(i)  <  %)  <  •  <  ^(n)> 

the  empirical  cumulative  distribution  function  is  [ 

{0,  x  <  x(i) 

r/n,  x{r)  <x<  x(r+i)  .  (9) 

1?  n )  —  % 

If  Fo(x)  is  the  true,  fully  specified  theoretical  cumulative  distribution  function  from  which  the  samples 
are  drawn,  then  from  the  strong  law  of  large  numbers 


lim  P{Sn(x)  =  F0(x)}  =  1. 

n — >oc 

Define  the  metric  used  for  measuring  the  separation  between  theoretical  and  empirical  CDFs  to  be 


Ar 


(F0(x)  -  Sn(x))dx 


This  is  different  from  the  usual  measure  of  deviation 

Dn  =  sup  |5n(x)  -  F0(x) | 


(10) 


(11) 


and  its  variants  that  are  well  developed  for  the  KS  test.  The  measure  of  deviation  (11)  proved  to  not 
be  sensitive  enough  for  the  data  association  problem  considered  here.  The  metric  (10)  is  more  sensitive  to 
incorrect  measurements.  This  results  from  the  fact  that  (10)  largely  takes  advantage  of  how  a  single  incorrect 
measurement  affects  the  shape  of  the  empirical  CDF.  When  using  ( 10)  distinct  peaks  appear  in  this  value 
for  windows  of  n  measurements  that  contain  false  measurements.  Figure  1  shows  these  peaks  in  error  for 
the  two  satellite  system  for  which  results  are  presented  in  Section  6. 


An 

jatHj 


Measurements 


200 


Figure  1:  Peaks  in  error  between  empirical  and  theoretical  CDFs  around  incorrect  measurements 
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Figure  2:  Satellites  in  close  proximity  with  measurement  mixing  occurring 


There  has  been  a  large  body  of  theory  developed  using  Dn  in  (1 1),  demonstrating  many  properties,  such 
as  the  distribution  of  Dn  is  independent  of  the  distribution  Fq(x)  for  continuous  CDFs  [  i].  Critical  values 
or  thresholds  dn(a)  have  also  been  established  [ .,  ]  such  that 

P{Dn  <  d„(a)}  =  1  -  a 

where  1  —  a  is  the  confidence  level.  That  is,  if  Dn  <  dn(a),  then  with  probability  1  —  a,  the  empirical 
CDF  Sn(x)  is  formed  from  n  samples  drawn  from  Fq(x).  However,  there  is  currently  no  theoretical  basis 
for  specifying  critical  values  for  An  in  (10),  and  this  is  an  area  of  future  work.  Section  (i  will  show  the 
performance  of  using  (11)  on  the  example  discussed  in  Section  5. 

4  Computational  Complexity 

The  batch  filter  has  historically  found  wide  use  throughout  the  aerospace  tracking  community.  A  probable 
reason  for  this  is  that  the  batch  fiber  converges  to  zero  RMS  error  with  the  assumption  of  a  perfect  model. 
The  batch  filter  has  been  used  classically  for  problems  where  there  is  no  measurement  mixing.  The  capability 
of  associating  data  is  not  directly  built  into  the  batch  filter  like  it  is  in  gating,  the  KS  tests,  and  PDA.  The 
batch  filter  can  be  used  for  the  data  association  problem,  but  at  a  large  cost  in  computational  complexity. 

To  illustrate  the  heavy  computational  complexity  of  applying  the  batch  filter  to  the  data  association 
problem,  consider  the  situation  pictured  in  Figure  2.  There  are  two  satellites  in  similar  orbits,  i.e.,  they  are 
in  close  proximity.  A  single  ground  based  sensor  produces  three  measurements  after  three  time  steps.  The 
first  and  third  measurements  come  from  satellite  A  and  the  second  from  satellite  B)  which  is  assumed  to 
not  be  of  interest.  The  method  by  which  the  batch  filter  accomplishes  data  association  is  to  minimize  the 
covariance  associated  with  the  estimate  over  the  set  of  possible  hypotheses  of  measurement  origin.  The  set 
of  possible  hypotheses  for  the  two  satellite  system  in  Figure  2  is 


7 


Hypothesis 

Measurement 

1  2  3 

a 

Y 

Y 

Y 

b 

Y 

Y 

N 

c 

Y 

N 

Y 

d 

Y 

N 

N 

e 

N 

Y 

Y 

f 

N 

Y 

N 

g 

N 

N 

Y 

h 

N 

N 

N 

Thus,  the  batch  filter  associates  the  three  measurements  by  testing  each  of  these  hypotheses  (i.e  running 
he  measurements  through  a  batch  filter  with  the  corresponding  assumptions)  to  find  the  one  that  minimizes 
Mnanance  o  he  estimate.  Obviously  for  this  case  the  covariance  will  be  minimized  by  hypothesis  c. 

T1  Fh  i  ,  mflmP  e  abv°Ve’  there  are  °n  y  tW°  satellites  and  three  measurements  taken  over  three  time  steps 
The  batch  filter  must  be  run  eight  separate  times  to  check  the  hypotheses  of  whether  or  not  each  individual 

a  SrmtTr  aTSch  d  T  ^  ^  °f  int^rest  "ot.  In  general,  the  computational  complexity  of  using 

a  batch  filter  at  each  individual  time  step  is  2fc,  where  k  denotes  the  current  number  of  time  steps  "that  have 

elapsed  since  measurements  started  to  be  taken.  For  the  example  above,  k  =  3.  It  is  easy  to  see  that  as  the 

thebatch^fflter  for  data  a  henCe,meafurements  Soes  UP> the  computational  complexity  associated  with  using 
tne  batch  hlter  tor  data  association  becomes  intractable  very  quickly. 

_  th,e  c°mPutational  complexity  of  hypothesis  testing  does  not  show  up  in  recursive 

Tods  ike  the  EKF  and  PDAF.  These  methods  only  take  the  current  measurement  as  input,  independent 
time  elapsed  since  starting  to  take  measurements  (i.e.,  mk  =  1,  the  number  of  measurements  used  for 

filtering  is  always  one  at  each  time  step).  Thus  there  is  a  constant  amount  of  computational  complexity 
associated  with  these  methods  at  each  time  step.  complexity 


5  Orbital  Description  *- 

To  actually  run  simulations  for  the  purposes  of  testing  these  data  association  techniques  orbits  for  the 

^  '  O  b/  “  ■ The  "bitS  *""•*  “  *•»•*  »'  th.  six  classical  dal® 

(  ,e.l,n ,u,  u)  denoting  semi-major  axis,  eccentricity,  inclination,  right  ascension  of  the  ascending  node 
argnment  o  perigee,  and  tine  anomaly.  The  snm  of  results  for  this  report  will  be  given  for  too  saSitet 

0nbin  a  tme  drCUl!r  °rbit  ^  the  °ther  WHh  a  Slight  P^turbation  in  the  eccentricity. 
Choosing  a  Cartesian  basis  corresponding  to  the  ECI  frame  provides  us  with  the  most  natural  choice 

of  coordinates  for  describing  orbits.  In  simulation,  the  six  orbital  parameters  are  used  to  obtain  position 
ve  ocity  initial  conditions  in  Cartesian  coordinates  aligned  with  the  ECI  frame.  Thus,  the  dynamics 
describing  the  orbits  are  written  down  in  this  frame.  The  method  by  which  the  dynamics  are  solved  is  a 
straightforward  approach  using  Lagrangian  techniques.  The  resulting  equations  are  simple  2-body  dynamics- 
the  same  equations  can  be  found  in  any  standard  astro-dynamics  text,  for  example  fi  1 

tlK  d5’’an,iCS  *"  ^  “5i"g  E^9ran9e, ovulation,.  The 

L  =  KB  —  PE.  (12) 

IZtTpT-  Afhe  ECI  !,rame’  t/-2killfC  energy  is  just  KE  =  1/2  Msatellitei  (i2  +  if  +  ?)  and  the  potential 
energy  rtj  —  Msateuite.Gearthy/x 2  +  y2  +  z2. 

The  dynamics  can  then  be  directly  obtained  from  the  Euler-Lagrange  equations 


where  q  is  the  state  vector  q  =  (x,x,p,p,  i,  2).  The  resulting  equations  arc 


x(t) 
x(t) 
v(t ) 
y(t) 

z(t) 

z{t) 


[ix(t) 


(x(t)2+y(t)2  +  z(t)2)^ 
m(t)  


(1  y  . 

sI(t) 


(x(t)2  +  y(t)i  +  z(t)*)y* 

d  ,  . 

dix(t) 

_ t^jt) _ 

(x(t)2  +  y{t)2  +  2(f)2)3/2 

d  y  . 

sI(t) 


where  p  =  3.986  •  105  kmV2.  These  equations  are  integrated  in  simulation  to  produce  the  “truth”  orbits, 
i.e.,  the  actual  orbit  with  initial  conditions  specified  by  the  six  orbital  parameters.  At  this  point  there  is  no 
noise  added  yet. 

The  sensor  is  located  on  the  surface  of  the  earth,  and  the  SEZ  frame  whose  origin  is  aligned  with  the 
sensor’s  location  is  the  frame  in  which  measurements  are  assumed  to  be  “taken.”  To  generate  “measurements” 
for  the  simulation,  the  positions  of  the  orbits  in  ECI  are  first  sampled  at  a  specified  rate.  Process  noise  is 
added  to  these  samples  because  the  dynamics  are  written  dowm  in  the  ECI  frame  and  the  process  noise 
is  the  noise  associated  with  inadequacies  in  this  model  (see  equation  (1)).  After  adding  process  noise,  a 
representation  of  the  resulting  “noisy”  ECI  samples  is  obtained  in  SEZ  coordinates  centered  at  the  sensor’s 
location  through  the  use  of  rigid  body  transformations. 

Obtaining  the  SEZ  representations  of  the  orbits  (with  process  noise  added)  is  accomplished  using  a 
composition  of  rigid  body  transformations.  The  resulting  transformation  has  the  form 


93  =  9 31  *  9i 

where  q\  is  the  homogeneous  representation  of  the  position  vector  q\  =  (x,  y ,  2, 1)  in  frame  1,  is  standard 
matrix  multiplication,  and  g  —  (p,  R)  £  SE( 3).  In  general,  g  can  be  written  as 


and  is  such  that 

£31  =  £32  •  £21  (14) 

where  R  is  a  standard  Euler  rotation  matrix  about  the  x,  p,  or  2  axes,  p  is  the  position  vector  (x,p,  2), 
and  is  again  standard  matrix  multiplication.  The  notation  g$i  defines  a  rigid  body  transform  g  that 
allows  us  to  transform  objects  represented  in  frame  1  to  frame  3.  The  right  hand  side  of  equation  (14)  is  a 
composition  of  rigid  body  transformations:  the  resulting  transformation  from  frame  1  to  frame  3  is  just  the 
transformation  from  frame  1  to  2  followed  by  the  transformation  from  frame  2  to  3.  For  a  more  thorough 
description  of  rigid  body  transformations  see 

To  transition  between  ECI  and  SEZ  coordinates  it  is  only  necessary  to  know  the  earth’s  rotational  rate, 
the  latitude  and  longitude  of  the  sensor’s  position  on  the  earth,  and  the  radius  of  the  earth.  These  four 
quantities  will  describe  four  rigid  body  transformations.  The  first  transformation  is  a  time  varying  rotation 
about  the  earth’s  rotation  axis  (which  in  ECI  is  aligned  with  the  2-axis)  at  a  rate  equal  to  the  earth’s 
rotation  rate,  and  aligns  the  frame  such  that  it  is  not  rotating  with  respect  to  the  sensor.  The  second 
transformation  is  another  rotation  about  the  earth’s  rotational  axis,  which  is  constant  and  corresponds  to 
the  sensor’s  longitude.  The  next  transformation  is  another  rotation,  this  one  about  the  p-axis  of  the  frame 
by  the  correct  amount  of  latitude  corresponding  to  the  sensor  position.  Finally,  the  last  transformation  is 
just  a  translation  in  the  2-direction  by  the  earth’s  radius.  Thus  the  resulting  rigid  body  transformation  that 
takes  points  represented  in  ECI  coordinates  into  SEZ  coordinates  is 
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9ECi^SEZ=9{Per,Rm-9(O,RyW-9(O,Rz(4>))-9(O,Rz(0m  (15) 

where  per  =  (0,0,6378.135  km),  Ry  an  Euler  rotation  matrix  corresponding  to  a  rotation  about  the  y-axis, 
'ip  the  latitude  at  which  the  sensor  is  located,  Rz  an  Euler  rotation  matrix  about  the-2-axis,  (p  the  longitude 
at  which  the  sensor  is  located,  and  0(t)  the  earth’s  rotational  rate. 

All  of  the  samples  produced  in  ECI  coordinates  undergo  this  transformation  and  we  are  left  with  (a;,  y,  z) 
descriptions  of  the  orbits  in  the  SEZ  frame.  These  transformations  can  clearly  be  seen  in  the  short  movie 
located  at  http://robotics.colorado.edU/wiki/index.php/Travers:Satellites. 

The  next  step  in  obtaining  measurements  for  the  simulation  is  to  then  take  the  noisy  samples  expressed 
in  an  SEZ  coordinate  frame  and  represent  them  in  R/Az/El  coordinates.  This  transformation  is  defined  by 

( R,Az,El )  =  {yj x2  +  y2  +  z2,  tan _1(— x/y),  sin”1  (z/\/x2  -fy2  +  z2)).  (16) 

Because  the  actual  sensors  take  measurements  in  R/Az/El,  measurement  noise  distributions  are  assumed  to 
be  represented  in  R/Az/El.  This  fact  implies  that  measurement  noise  is  added  to  the  samples  in  R/Az/El. 

These  measurements  in  R/Az/El  represented  in  a  sensor  local  frame  then  need  to  be  transformed  back 
into  ECI  coordinates  so  that  they  are  in  the  form  of  equation  (2).  This  is  accomplished  by  doing  the 
previously  stated  transformations  in  reverse  order.  Note  that  the  inverse  transformation  for  equation  (16)  is 

(a;,  y,  z)  =  (— Rcos(Az)  cos  (El),  Rsin(Az)  cos(£7),  Rsin(£7)),  (17) 

and  the  rigid  body  transformation  from  SEZ  coordinates  to  ECI  coordinates  is  the  matrix  product  on  the 
right  hand  side  of  (15)  in  reverse  order. 

The  fact  that  measurements  are  taken  in  R/Az/El  coordinates  from  an  SEZ  coordinate  frame  implies 
that  the  resulting  measurement  error  covariance  is  most  naturally  represented  in  R/Az/El,  i.e.,  the  measure¬ 
ment  noise  distributions  in  all  three  directions  R/Az/El  are  Gaussian.  The  measurement  error  covariance, 
RR/Az/Ei(k)i  must  be  represented  in  the  ECI  frame  to  correspond  with  the  measurement  and  system  model. 
The  variance  in  each  direction  can  be  thought  of  as  a  member  of  the  corresponding  tangent  space  at  each 
point,  thus  the  transformation  of  the  covariance  has  the  form 

(fj(*>)  •  (W(‘)l  ■  (§j(*>)' ’  (18) 

where  dx/dy  is  the  Jacobian.  In  R/Az/El  coordinates  the  measurement  error  covariance  description  is 
constant,  under  this  transformation  the  measurement  error  covariance  becomes  time  varying.  A  movie  of  how 
the  three  dimensional  measurement  error  distribution  changes  as  an  object  traverses  an  orbit,  when  viewed 
from  ECI  coordinates,  can  be  found  at  http://rohotics.colorado.edU/wiki/index.php/Travers:Satellites. 


6  Results 

Results  are  presented  for  a  simulated  system  consisting  of  two  satellites.  A  number  of  system  parameters  are 
varied  to  highlight  performance  over  two  specific  objectives:  tracking  capability  and  correct  identification  of 
individual  measurements  as  being  from  the  satellite  being  tracked  or  not. 

Figure  3  displays  the  results  of  each  of  the  filters  described  in  Section  3  in  terms  of  RMS  position  tracking 
error.  The  RMS  error  is  the  error  between  the  filtered  estimate  of  position  at  each  discrete  time  and  the 
actual  position  of  the  simulated  satellite  in  orbit  at  each  corresponding  time.  The  cumulative  average  of 
the  norms  (standard  Euclidean  norm)  of  theses  differences  up  to  each  time  step  are  the  values  of  the  RMS 
calculation. 

The  EKF  without  false  measurements  as  well  as  the  Batch  filter  are  run  on  a  list  of  measurements  that  do 
not  include  mixing  (all  the  measurements  are  from  the  same  satellite).  These  results  are  included  as  a  basis 
for  comparison.  The  same  can  be  said  about  the  “Samples  Alone”  plot.  This  plot  is  obtained  by  comparing 
the  value  of  the  position  measurement  of  a  single  satellite  expressed  in  ECI  coordinates  (with  process  noise 
and  measurement  noise  added  in  the  appropriate  frames)  to  truth  at  each  time  step,  then  calculating  the 
new  RMS  value  including  the  cumulative  normed  difference  up  to  that  time  step.  This  plot  can  be  seen  as 
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Figure  3:  Comparison  of  various  filter  performances.  The  horizontal  axis  is  time  in  seconds  and  the  vertical 
RMS  error  in  kilometers.  The  upper  plot  presents  results  for  the  case  where  the  two  satellites  are  further 
away  (e/  =  0.01)  from  each  other  than  the  lower  (ec  =  0.0004).  The  UI=95%. 


how  well  the  measurements  from  a  single  satellite  do  without  any  filtering.  In  terms  of  RMS  error  of  tracking 
a  single  satellite  where  there  is  no  measurement  mixing,  the  batch  filter  performs  better  than  the  EKF.  Both 
methods  offer  an  improvement  in  tracking,  measured  in  RMS  error,  over  the  measurements  considered  alone. 

The  single  satellite  for  which  these  results  are  stated  is  in  a  circular  equatorial  orbit.  The  radius  of 
this  orbit  is  r  =  7078.14  km.  The  initial  conditions  are  always  such  that  the  satellite  starts  on  the  a:- axis 
of  the  ECI  frame.  The  results  of  the  other  filters  evaluated  here  are  all  for  a  two  satellite  system.  The 
uniformity  index  is  95%,  meaning  that  5%  of  the  total  measurements  are  from  satellite  2,  while  it  is  satellite 
1  that  we  are  trying  to  track.  Proximity  for  these  two  satellites  will  be  defined  in  terms  of  a  perturbation 
in  eccentricity.  The  first  satellite  is  in  the  same  orbit  described  above,  i.e. ,  a  circular  equatorial  orbit  with 
r  =  7078.14  km.  The  second  satellite  is  again  equatorial,  except  that  it  is  slightly  elliptical  with  semi-major 
axis  a  =  7078.14  km.  For  the  “Far-”  proximity  case,  ej  =  0.01  and  for  the  “Close”  case  ec  —  0.0004.  Over  the 
measurement  period  of  2s,  Ef  corresponds  to  an  average  separation  distance  of  70.78  km  and  ec  corresponds 
to  an  average  separation  of  8.9  km.  The  process  noise  associated  with  each  component  of  the  position  is 
described  by  the  normal  distribution  N( 0,  (0.7)2km2).  The  process  noise  associated  with  each  component  of 
the  velocity  is  described  by  the  distribution  1V( 0,  (3)2(km/s)2).  The  distributions  describing  measurement 
noise  in  range,  azimuth,  and  elevation  are  N( 0,  (0.02)2km2),  N( 0,  (0.01)2radians2),  and  JV(0,  (0.01)2radians2) 
correspondingly.  Lastly,  the  time  step  is  dt  =  0.001  s. 

The  PDA  results  are  in  general  better  than  the  EKF  and  converge  faster  than  the  Batch  filter.  The 
EKF  with  Gating  does  better  than  the  EKF  for  the  no  mixing  case  as  well  as  the  EKF  run  on  the  results  of 
the  KS-Test.  When  the  two  satellites  are  far  from  each  other,  the  EKF  with  Gating  does  noticeably  better 
than  the  other  two  EKF  results,  while  the  EKF  with  KS  does  basically  about  as  good  as  the  EKF  evaluated 
without  any  incorrect  measurements.  The  EKF  with  Gating  results  do  better  because  at  this  distance  they 
not  only  gate  out  all  of  the  incorrect  measurements,  but  they  also  gate  out  the  true  measurements  that  are 
in  effect  bad  measurements.”  The  EKF  run  on  KS  results  and  the  EKF  on  true  measurements  alone  are 
essentially  the  same  because  the  KS  test  only  removes  the  incorrect  measurements,  meaning  that  after  the 
KS  test  is  run,  the  resulting  list  of  measurements  is  exactly  the  same  as  the  list  of  measurements  used  to 
evaluate  the  EKF  alone,  except  that  in  this  case  5%  of  the  measurements  are  removed  by  the  test. 

When  the  two  satellites  are  in  close  proximity  the  EKF  results  change  slightly.  The  EKF  run  on  true 
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Figure  4:  Comparison  of  PDA  with  gating  performance  as  UI  is  varied.  RMS  error  in  kilometers  are  plotted 
against  time  in  seconds.  The  satellites  considered  are  both  in  equatorial  orbits  with  a  =  7078.14.  The  satellite 
being  tracked  has  eccentricity,  £t  =  0.1  and  the  satellite  not  being  tracked  has  eccentricity,  £fi  =0.11. 


measurements  obviously  docs  not  change.  Both  the  EKF  with  Gating  and  EKF  with  KS  test  do  worse  as 
the  satellites  come  into  closer  proximity.  The  EKF  with  gating  still  does  marginally  better  than  the  EKF 
on  truth,  but  the  KS  test  results  are  actually  worse.  This  can  be  explained  by  looking  at  Figure  5.  At 
the  separation  distance  of  8.9  km,  the  KS  test  is  almost  completely  failing.  Thus  most  of  the  incorrect 
measurements  are  being  left  in  the  list  of  measurements  and  the  filtering  performance  is  adversely  affected. 

Figure  4  illustrates  the  effect  of  varying  the  uniformity  index  (UI)  on  PDA  filtering  performance.  As 
the  UI  goes  down,  filtering  performance  degrades.  For  the  case  plotted,  two  satellites  in  equatorial  orbits 
are  considered  with  a  =  7078.14  km  for  each.  The  satellite  being  tracked  has  eccentricity,  et  =  0.1,  and  the 
“follower”  satellite  has  eccentricity,  eji  =0.11. 

Described  in  Section  4,  both  gating  and  the  KS  test  variant  are  methods  of  identifying  individual  mea¬ 
surements  as  being  from  the  satellite  being  tracked  or  not.  Figure  5  displays  simulated  results  for  the  two 
satellite  system  where  the  separation  distance  between  the  two  satellites  is  varied  and  the  corresponding 
percentage  of  incorrect  measurements  correctly  identified  is  calculated  for  each  method.  The  separation 
distance  is  calculated  from  the  average  separation  between  satellites  over  the  total  measurement  period.  In 
simulation  the  separation  distance  is  actually  varied  by  varying  the  eccentricity  of  the  orbit  of  satellite  2, 
other  than  this  the  system  parameters  are  the  same  as  above  with  UI  =  95%.  The  KS  test  starts  to  have 
degraded  performance  at  a  higher  separation  distance  than  docs  gating,  although  above  about  40  km  both 
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Figure  5:  Comparison  of  Gating  and  KS  methods  at  correctly  identifying  false  measurements  as  the  separa¬ 
tion  distance  between  the  two  satellites  is  varied;  UI=95%. 
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Figure  6:  Comparison  of  Gating  and  KS  variant  as  a  function  of  the  discrepancy  between  actual  and  modeled 
measurement  error  covariances.  The  UI=95%.  The  two  satellites  considered  are  in  the  “far”  configuration 
(i.e.,  £/  =  0.01). 


methods  perfectly  identify  all  of  the  false  measurements. 

Correctly  identifying  incorrect  measurements  carries  a  tradeoff  in  gating:  incorrectly  identifying  true 
measurements  as  false.  The  intent  of  this  work  is  not  only  to  correctly  identify  false  measurements  as  false, 
but  also  to  correctly  identify  true  measurements  as  true  while  doing  so.  Figure  0  shows  a  particular  situation 
in  which  gating  struggles  to  correctly  identify  true  measurements.  For  the  case  considered,  UI  =  95%  and  the 
two  satellites  are  in  the  far  configuration  (i.e.,  the  two  satellites  are  in  equatorial  orbits,  the  satellite  being 
tracked  in  a  circular  orbit  and  the  satellite  not  being  tracked  having  identical  orbital  element  values  except 
for  a  difference  in  eccentricity,  £f  =  0.01).  This  figure  plots  the  percentage  of  true  measurements  correctly 
identified  as  true  as  a  function  of  how  well  the  modeled  and  actual  measurement  error  covariances  match 
each  other.  Thus,  at  Sigma=l  the  modeled  and  actual  covariance’s  match  each  other  perfectly.  At  Sigma=2 
the  actual  measurement  error  covariance  is  twice  the  modeled,  and  so  on.  As  the  discrepancy  between  actual 
and  modeled  measurement  error  covariance  increases,  the  performance  of  gating  can  be  seen  to  degrade 
quickly,  while  the  KS  test  maintains  its  performance  throughout.  This  is  explained  by  the  empirical  basis 
from  which  the  KS  test  can  form  “theoretical”  CDFs.  The  KS  test  does  not  rely  on  an  accurate  model  of 
the  measurement  error  covariance.  Rather,  when  the  satellite  of  interest  is  not  close  to  other  space  objects, 
“clean”  unmixed  measurements  can  be  used  to  empirically  determine  the  “theoretical”  CDF  for  the  satellite. 
Hence,  the  KS  test  also  does  not  depend  on  Gaussian  distributions  of  the  measurement  errors. 


7  Conclusions  &:  Future  Works 

This  report  has  highlighted  several  key  issues  facing  tracking  and  correct  data  association  for  multiple  satellite 
systems  where  there  is  measurement  mixing.  Initial  solutions  to  these  issues  were  presented,  and  simulation 
results  given. 

The  first  problem  addressed  was  that  of  computational  complexity,  a  common  concern  in  typical  tracking 
schemes  for  multi-body  systems.  The  historical  tendency  to  use  batch  filters  for  the  purposes  of  filtering 
single  body  systems  turns  out  to  be  computationally  burdensome  for  multi-body  systems  with  measurement 
mixing.  As  stated,  the  batch  approach  to  data  association  grows  as  2k  where  k  is  the  number  of  time 
steps.  For  one  measurement  (of  unknown  origin)  at  each  time  step  while  tracking  one  satellite,  the  recursive 
methods  such  as  the  EKF  and  PDAF  outlined  in  Section  5  have  constant  complexity.  For  this  reason  the 
recursive  filtering  methods  are  preferable  for  filtering  satellite  systems  where  there  is  measurement  mixing. 

Simulation  results  where  presented  for  a  simple  two  satellite  system  with  measurement  mixing.  The 
results  where  grouped  into  filtering  performance  and  correct  data  association  performance.  The  results  show 
that  of  the  data  association  and  filtering  algorithms  considered,  PDA  with  gating  has  the  best  performance 
although  PDA  by  itself  still  offers  a  dramatic  improvement  over  the  standard  EKF  used  by  itself  or  in 
conjunction  with  the  gating  or  KS  techniques. 
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The  data  association  problem  was  considered  from  two  separate  standpoints,  the  first  being  the  identifica¬ 
tion  of  false  measurements  and  the  second  the  correct  identification  of  true  measurements  while  identifying 
false  measurements.  In  terms  of  finding  a  minimum  separation  distance  between  satellites  for  which  the 
identification  of  incorrect  measurements  is  possible,  gating  is  the  preferred  method  of  data  association  over 
the  KS  test  variant  In  terms  of  correctly  identifying  true  measurements,  the  performance  of  gating  contains 
tradeoffs  and  in  general  is  not  very  robust  to  modeling  errors  in  the  system  parameters.  This  was  highlighted 
in  the  simulation  results  where  the  difference  between  modeled  and  actual  measurement  error  covariances  for 
the  simulated  system  was  varied.  The  empirical  basis  of  the  KS  test  variant  allowed  for  perfect  performance 
as  the  modeled  covariance  deviated  further  from  the  truth,  while  the  gating  method  performance  degraded. 
These  results  show  that  both  gating  and  the  KS  test  variant  have  favorable  characteristics  in  different  situ¬ 
ations.  For  this  reason,  one  proposed  avenue  of  future  work  is  the  development  of  a  hybrid  method  of  these 
two  techniques  such  that  the  strongest  qualities  of  each  are  used  to  produce  a  better  data  association  method 
for  systems  with  mixed  measurements. 

Another  focus  of  future  work  is  to  add  various  levels  of  complexity  and  fidelity  to  the  two  satellite, 
single  sensor  system  considered  in  the  Section  0.  The  first  goal  is  to  implement  a  network  of  sensors  around 
the  earth  that  take  measurements  of  the  satellites  as  they  travel  into  and  out  of  each  individual  sensor’s 
range  (which  is  limited  to  line  of  sight).  The  work  presented  in  this  report  assumes  that  line  of  sight  is 
not  required  for  the  sensor  to  make  measurements  (i.e.,  the  single  sensor  on  the  surface  of  the  earth  is 
assumed  to  be  able  to  track  a  satellite  through  an  entire  orbit).  The  network  of  sensors  is  a  more  realistic 
scheme  and  in  theory  could  present  added  complications  to  the  data  association  problem,  especially  in  regions 
where  the  object  being  tracked  switches  from  one  sensor’s  range  to  another.  Adding  maneuvering  to  the 
individual  satellites  is  another  area  of  future  work.  With  maneuvering,  the  data  association  problem  can  be 
studied  for  systems  where  two  satellites  are  originally  not  in  close  proximity,  but  then  one  of  the  satellites 
moves  toward  the  other  (“chasing”  it  in  a  sense).  In  addition  to  a  network  of  sensors  and  maneuvering, 
the  effect  of  adding  more  satellites  to  the  system  will  be  considered  in  future  work.  The  results  presented 
in  this  report  considered  two  satellites  with  proximity  defined  in  a  specific  way  (namely  the  two  by  one 
ellipse  that  results  from  a  perturbation  in  eccentricity  of  two  satellites  in  the  same  orbit  with  the  same  initial 
conditions).  The  techniques  presented  in  this  report  are  extendable  to  multi-satellite  systems,  but  the  overall 
level  of  complexity  that  results  from  having  many  satellites  in  varied  degrees  of  proximity  may  require  more 
advanced  data  association  techniques. 
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